#define _H_Mesh_Private

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "Mesh.h"
#include "PMesh.h"
#include "MSTK.h"
#include "MSTK_private.h"
#include "MSTK_malloc.h"


#ifdef __cplusplus
extern "C" {
#endif

static  RepType MESH_rtype[5] = {F1, F4, R1, R2, R4};
static  char MESH_rtype_str[5][3] = {"F1\0","F4\0","R1\0","R2\0","R4\0"};

  /*for mpi*/
  static  char MESH_ptype_str[4][12] = {"PINTERIOR\0","PUTILIZED\0","PBOUNDARY\0","PGHOST\0"};
  /*end for mpi*/

  /* this function convert a regular mstk vertex into PVertex type for data transfer */


  /* this function copies information from mv to mv_new 
   mv_new should be allocated before this call:
   mv_new = MV_New(mesh)
   PType is not copied
  */
void PMV_Duplicate(MVertex_ptr mv, MVertex_ptr mv_new) {
  double  xyz[3];
  MV_Coords(mv,xyz);
  MV_Set_Coords(mv_new,xyz);
  MV_Set_GEntID(mv_new,MV_GEntID(mv));
  MV_Set_GEntDim(mv_new,MV_GEntDim(mv));
  MV_Set_MasterParID(mv_new,MV_MasterParID(mv));
  MV_Set_GlobalID(mv_new,MV_GlobalID(mv));
}

void PME_Duplicate(MEdge_ptr me, MEdge_ptr me_new) {
  ME_Set_GEntID(me_new,ME_GEntID(me));
  ME_Set_GEntDim(me_new,ME_GEntDim(me));
  ME_Set_MasterParID(me_new,ME_MasterParID(me));
  ME_Set_GlobalID(me_new,ME_GlobalID(me));
}

void PMF_Duplicate(MFace_ptr mf, MFace_ptr mf_new) {
  MF_Set_GEntID(mf_new,MF_GEntID(mf));
  MF_Set_GEntDim(mf_new,MF_GEntDim(mf));
  MF_Set_MasterParID(mf_new,MF_MasterParID(mf));
  MF_Set_GlobalID(mf_new,MF_GlobalID(mf));
}


void PMR_Duplicate(MRegion_ptr mr, MRegion_ptr mr_new) {
  MR_Set_GEntID(mr_new,MR_GEntID(mr));
  MR_Set_GEntDim(mr_new,MR_GEntDim(mr));
  MR_Set_MasterParID(mr_new,MR_MasterParID(mr));
  MR_Set_GlobalID(mr_new,MR_GlobalID(mr));
}



PType Str2PType(char *str) {
  if (strncmp(str,"PINTERIOR",3) == 0)
      return PINTERIOR;
  if (strncmp(str,"PUTILIZED",3) == 0)
      return PUTILIZED;
  if (strncmp(str,"PBOUNDARY",3) == 0)
      return PBOUNDARY;
  if (strncmp(str,"PGHOST",3) == 0)
      return PGHOST;
  MSTK_Report("Str2PType","PType unmatch",FATAL);
  return 0;
}

  /* File name should have the .mstk extension - will not check here */
int PMESH_InitFromFile(Mesh_ptr mesh, const char *filename) {
  FILE *fp;
  char inp_rtype[16], temp_str[256], fltype_str[16], rltype_str[16];
  char attname[256], atttype_str[256], attent_str[256];
  int i, j, found, NV=0, NE=0, NF=0, NR=0, nav, nar, gdim, gid;
  int dim, vid1, vid2, eid, fid, rid, adjvid, adjrid, adjv_flag;
  int nfv, max_nfv=0, nfe, max_nfe=0, nrv, max_nrv=0, nrf, max_nrf=0;
  int *fedirs, *rfdirs, ival, ncomp, nent, done, status;
  int processed_vertices=0, processed_adjv=0, processed_edges=0;
  int processed_faces=0, processed_regions=0, processed_adjr=0;
  double ver, xyz[3], rval;
  double *rval_arr;
  MVertex_ptr mv, ev1, ev2, adjv, *fverts, *rverts;
  MEdge_ptr me, *fedges;
  MFace_ptr mf, *rfaces;
  MRegion_ptr mr, adjr;
  MEntity_ptr ent;
  MAttrib_ptr attrib;
  MType attent;
  MAttType atttype;
  /* for mpi */
  char ptype_str[10];
  int id, masterparid, globalid;
  /* end for mpi */

  if (!(fp = fopen(filename,"r"))) {
    MSTK_Report("MESH_InitFromFile","Cannot open file",ERROR);
    return 0;
  }

  status = fscanf(fp,"%s %lf",temp_str,&ver);
  if (strcmp(temp_str,"MSTK") != 0) {
    MSTK_Report("MESH_InitFromFile","Not a MSTK file",ERROR);
    fclose(fp);
    return 0;
  }
  if (status == EOF)
    MSTK_Report("MESH_InitFromFile",
		"Premature end of file before any mesh data is read",FATAL);

  if (ver != MSTK_VER) {
    MSTK_Report("MESH_InitFromFile","Version mismatch",WARN);
  }

  status = fscanf(fp,"%s %d %d %d %d\n",inp_rtype,
		  &(mesh->nv),&(mesh->ne),&(mesh->nf),&(mesh->nr));
  if (status == EOF)
    MSTK_Report("MESH_InitFromFile",
		"Premature end of file before any mesh data is read",FATAL);


  found = 0;
  for (i = 0; i < MSTK_MAXREP; i++) {
    if (strncmp(inp_rtype,MESH_rtype_str[i],2) == 0) {
      MESH_SetRepTypeIni(mesh, MESH_rtype[i]);
      found = 1;
      break;
    }
  }

  if (!found) {
    MSTK_Report("MESH_InitFromFile","Unrecognized representation type",ERROR);
    fclose(fp);
    return 0;
  }


  /* For now, the reduced representations only allow region and
     face elements - no edge elements are allowed */

  if (mesh->reptype >= R1 && mesh->reptype <= R4) {
    if (mesh->ne)
      MSTK_Report("Mesh_InitFromFile",
		  "Representation does not allow edges",WARN);
  }


  status = fscanf(fp,"%s",temp_str);
  if (status == EOF)
    MSTK_Report("MESH_InitFromFile",
		"Premature end of file while looking for vertex data",FATAL);
  else if (status == 0)
    MSTK_Report("MESH_InitFromFile",
		"Error in reading vertex data",FATAL);

  if (strncmp(temp_str,"vertices",8) == 0) {

    NV = mesh->nv;
    mesh->mvertex = List_New(NV);

    for (i = 0; i < NV; i++) {
      status = fscanf(fp,"%lf %lf %lf %d %d %s %d %d %d",&xyz[0],&xyz[1],&xyz[2],&gdim,&gid,
		      ptype_str,&id,&masterparid,&globalid);
      if (status == EOF)
	MSTK_Report("MESH_InitFromFile",
		    "Premature end of file while reading vertices",FATAL);
      else if (status == 0)
	MSTK_Report("MESH_InitFromFile",
		    "Error in reading vertex data",FATAL);

      mv = MV_New(mesh);
      MV_Set_Coords(mv,xyz);

      MV_Set_GEntID(mv,gid);      
      MV_Set_ID(mv,id);
      MV_Set_PType(mv,Str2PType(ptype_str));
      MV_Set_MasterParID(mv,masterparid);
      MV_Set_GlobalID(mv,globalid);
      MV_Set_GEntDim(mv,gdim);
      

    }

    processed_vertices = 1;
  }
  else {
    MSTK_Report("MESH_InitFromFile",
		"Vertex information should be listed first",ERROR);
    fclose(fp);
    return 0;
  }


  status = fscanf(fp,"%s",temp_str);
  if (status == EOF) {
    if (mesh->ne || mesh->nf || mesh->nr)
      MSTK_Report("MESH_InitFromFile",
		  "Premature end of file after vertex data",FATAL);
    else
      return 0;
  }



  /* ADJACENT VERTEX DATA */

  if (strncmp(temp_str,"adjvertices",11) == 0) {
    adjv_flag = 1;
    if (mesh->reptype == R2 || mesh->reptype == R4) {
      for (i = 0; i < NV; i++) {
	mv = List_Entry(mesh->mvertex, i);
	status = fscanf(fp,"%d",&nav);
	if (status == EOF)
	  MSTK_Report("MESH_InitFromFile",
		      "Premature end of file while reading adjacent vertices",
		      FATAL);
	else if (status == 0)
	  MSTK_Report("MESH_InitFromFile",
		      "Error in reading adjacent vertex data",FATAL);

	for (j = 0; j < nav; j++) {
	  status = fscanf(fp,"%d",&adjvid);
	  if (status == EOF)
	    MSTK_Report("MESH_InitFromFile",
			"Premature end of file while reading adjacent vertices"
			,FATAL);
	  else if (status == 0)
	    MSTK_Report("MESH_InitFromFile",
			"Error in reading adjacent vertex data",FATAL);

	  adjv = List_Entry(mesh->mvertex,adjvid-1);

	  
	  MV_Add_AdjVertex(mv,adjv);
	}
      }
    }
    else {
      for (i = 0; i < NV; i++) {
	status = fscanf(fp,"%d",&nav);
	if (status == EOF)
	  MSTK_Report("MESH_InitFromFile",
		      "Premature end of file while reading adjacent vertices",
		      FATAL);
	else if (status == 0)
	  MSTK_Report("MESH_InitFromFile",
		      "Error in reading adjacent vertex data",FATAL);
	
	for (j = 0; j < nav; j++)
	  fscanf(fp,"%d",&adjvid);
      }
    }

    processed_adjv = 1;
  }
  else {
    adjv_flag = 0;
    if (mesh->reptype == R2 || mesh->reptype == R4) {
      MSTK_Report("MESH_InitFromFile",
		  "Expected adjacent vertex information",ERROR);
    }
  }
    

  if (processed_adjv) {
    status = fscanf(fp,"%s",temp_str);
    if (status == EOF) {
      if (mesh->ne || mesh->nf || mesh->nr)
	MSTK_Report("MESH_InitFromFile",
		    "Premature end of file after adjacent vertex data",FATAL);
      else {
        fclose(fp);
	return 1;
      }
    }
  }



  /* EDGE DATA */

  if (strncmp(temp_str,"edges",5) == 0) {
 
    NE = mesh->ne;
    mesh->medge = List_New(NE);

    if (mesh->reptype >= F1 && mesh->reptype <= F4) {
      for (i = 0; i < NE; i++) {
	status = fscanf(fp,"%d %d %d %d %s %d %d %d",&vid1,&vid2,&gdim,&gid,
			ptype_str,&id,&masterparid,&globalid);
	if (status == EOF)
	  MSTK_Report("MESH_InitFromFile",
		      "Premature end of file while reading edges",FATAL);
	else if (status == 0)
	  MSTK_Report("MESH_InitFromFile",
		      "Error in reading edge data",FATAL);

	ev1 = List_Entry(mesh->mvertex,vid1-1);
	ev2 = List_Entry(mesh->mvertex,vid2-1);


	me = ME_New(mesh);

	ME_Set_Vertex(me,0,ev1);
	ME_Set_Vertex(me,1,ev2);


	ME_Set_GEntID(me,gid);
	ME_Set_ID(me, id);
	ME_Set_PType(me,Str2PType(ptype_str));
	ME_Set_MasterParID(me,masterparid);
	ME_Set_GlobalID(me,globalid);

	ME_Set_GEntDim(me,gdim);


      }
    }
    else {
      for (i = 0; i < NE; i++) {
	status = fscanf(fp,"%d %d %d %d %s %d %d %d",&vid1,&vid2,&gdim,&gid,
			ptype_str,&id,&masterparid,&globalid);
	if (status == EOF)
	  MSTK_Report("MESH_InitFromFile",
		      "Premature end of file while reading edges",FATAL);
	else if (status == 0)
	  MSTK_Report("MESH_InitFromFile",
		      "Error in reading edge data",FATAL);
      }
    }

    processed_edges = 1;
  }
  else {
    if (mesh->reptype >= F1 && mesh->reptype <= F4) {
      MSTK_Report("MESH_InitFromFile","Expected edge information",ERROR);
      return 0;
    }
  }
  
  if (processed_edges) {
    status = fscanf(fp,"%s",temp_str);
    if (status == EOF) {
      if (mesh->nf || mesh->nr) {
	MSTK_Report("MESH_InitFromFile",
		    "Premature end of file after edge data",FATAL);
      }
      else {
        fclose(fp);
	return 1;
      }
    }
  }



  /* FACE DATA */

  if (strncmp(temp_str,"face",4) == 0) {
    if (strncmp(temp_str,"faces",5) != 0) 
      MSTK_Report("MESH_InitFromFile","Expected keyword \"faces\"",ERROR);
    
    NF = mesh->nf;
    mesh->mface = List_New(NF);

    status = fscanf(fp,"%s",fltype_str);
    if (status == EOF)
      MSTK_Report("MESH_InitFromFile",
		  "Premature end of file while reading faces",FATAL);

    if (strncmp(fltype_str,"vertex",6) == 0) {      
      if (mesh->reptype >= R1 && mesh->reptype <= R4) {

	fverts = NULL;
	for (i = 0; i < NF; i++) {
	  mf = MF_New(mesh);

	  status = fscanf(fp,"%d",&nfv);
	  if (status == EOF)
	    MSTK_Report("MESH_InitFromFile",
			"Premature end of file while reading faces",FATAL);
	  else if (status == 0)
	    MSTK_Report("MESH_InitFromFile",
			"Error in reading face data",FATAL);

	  if (fverts) {
	    if (nfv > max_nfv) {
	      max_nfv = nfv;
	      fverts = MSTK_realloc(fverts,max_nfv*sizeof(MVertex_ptr));
	    }
	  }
	  else {
	    max_nfv = nfv;
	    fverts = MSTK_malloc(nfv*sizeof(MVertex_ptr));
	  }
	  
	  for (j = 0; j < nfv; j++) {
	    status = fscanf(fp,"%d",&vid1);
	    if (status == EOF)
	      MSTK_Report("MESH_InitFromFile",
			  "Premature end of file while reading face data",
			  FATAL);
	    else if (status == 0)
	      MSTK_Report("MESH_InitFromFile",
			  "Error in reading edge data",FATAL);

	    fverts[j] = List_Entry(mesh->mvertex,vid1-1);
	  }
	  MF_Set_Vertices(mf,nfv,fverts);  

	  status = fscanf(fp,"%d %d %s %d %d %d",&gdim,&gid,
			ptype_str,&id,&masterparid,&globalid);
	  if (status == EOF)
	    MSTK_Report("MESH_InitFromFile",
			"Premature end of file while reading face data",FATAL);
	  else if (status == 0)
	    MSTK_Report("MESH_InitFromFile",
			"Error in reading face data",FATAL);
	  
	  MF_Set_GEntID(mf,gid);
	  MF_Set_ID(mf,id);
	  MF_Set_PType(mf,Str2PType(ptype_str));
	  MF_Set_MasterParID(mf,masterparid);
	  MF_Set_GlobalID(mf,globalid);
	  MF_Set_GEntDim(mf,gdim);

	}
	if (fverts)
	  MSTK_free(fverts);

	processed_faces = 1;
      }
      else {
	MSTK_Report("MESH_InitFromFile",
		    "Expected face description in terms of vertices",ERROR);
        fclose(fp);
	return 0;
      }
    }
    else if (strncmp(fltype_str,"edge",4) == 0) {
      if (mesh->reptype >= F1 && mesh->reptype <= F4) { 

	fedges = NULL;	fedirs = NULL;
	for (i = 0; i < NF; i++) {
	  mf = MF_New(mesh);

	  status = fscanf(fp,"%d",&nfe);
	  if (status == EOF)
	    MSTK_Report("MESH_InitFromFile",
			"Premature end of file while reading face edge data",
			FATAL);
	  else if (status == 0)
	    MSTK_Report("MESH_InitFromFile",
			"Error in reading face data",FATAL);

	  if (fedges) {
	    if (nfe > max_nfe) {
	      max_nfe = nfe;
	      fedges = MSTK_realloc(fedges,max_nfe*sizeof(MVertex_ptr));
	      fedirs = MSTK_realloc(fedirs,max_nfe*sizeof(int));
	    }
	  }
	  else {
	    max_nfe = nfe;
	    fedges = MSTK_malloc(max_nfe*sizeof(MVertex_ptr));
	    fedirs = MSTK_malloc(max_nfe*sizeof(MVertex_ptr));
	  }
	  
	  for (j = 0; j < nfe; j++) {
	    status = fscanf(fp,"%d",&eid);
	    if (status == EOF)
	      MSTK_Report("MESH_InitFromFile",
			  "Premature end of file while reading face edge data",
			  FATAL);
	    else if (status == 0)
	      MSTK_Report("MESH_InitFromFile",
			  "Error in reading face data",FATAL);

	    fedirs[j] = eid > 0 ? 1 : 0;
	    fedges[j] = List_Entry(mesh->medge,abs(eid)-1);
	  }
	  
	  MF_Set_Edges(mf,nfe,fedges,fedirs);
	  
	  status = fscanf(fp,"%d %d %s %d %d %d",&gdim,&gid,
			  ptype_str,&id,&masterparid,&globalid);
	  if (status == EOF)
	    MSTK_Report("MESH_InitFromFile",
			"Premature end of file while reading faces",FATAL);
	  else if (status == 0)
	    MSTK_Report("MESH_InitFromFile",
			"Error in reading face data",FATAL);
	  
	  MF_Set_GEntDim(mf,gdim);
	  MF_Set_ID(mf,id);
	  MF_Set_PType(mf,Str2PType(ptype_str));
	  MF_Set_MasterParID(mf,masterparid);
	  MF_Set_GlobalID(mf,globalid);
	 
	  MF_Set_GEntID(mf,gid);
	}
	if (fedges) {
	  MSTK_free(fedges);
	  MSTK_free(fedirs);
	}
	
	processed_faces = 1;
      }
      else {
	MSTK_Report("MESH_InitFromFile",
		    "Expect face description in terms of vertices",ERROR);
	return 0;
      }
    }
  }
  else {
    if (mesh->reptype >= F1 && mesh->reptype <= F4) {
      MSTK_Report("MESH_InitFromFile","Expected face information",ERROR);
      fclose(fp);
      return 0;
    }
  }
  
  
  if (processed_faces) {
    status = fscanf(fp,"%s",temp_str);
    if (status == EOF) {
      if (mesh->nr != 0)
	MSTK_Report("MESH_InitFromFile",
		    "Premature end of file after face data",FATAL);
      else {
        fclose(fp);
	return 1;
      }
    }
    else if (status == 0)
      MSTK_Report("MESH_InitFromFile",
		  "Error in reading region data",FATAL);
  }


  /* REGION DATA */

  if (strncmp(temp_str,"region",6) == 0) {
    if (strncmp(temp_str,"regions",7) != 0) 
      MSTK_Report("MESH_InitFromFile","Expected keyword \"regions\"",ERROR);

    NR = mesh->nr;
    mesh->mregion = List_New(NR);

    status = fscanf(fp,"%s",rltype_str);
    if (status == EOF)
      MSTK_Report("MESH_InitFromFile",
		  "Premature end of file while reading regions",FATAL);
    else if (status == 0)
      MSTK_Report("MESH_InitFromFile",
		  "Error in reading region data",FATAL);

    if (strncmp(rltype_str,"vertex",6) == 0) {
      if (mesh->reptype == R1 || mesh->reptype == R2) {
	rverts = NULL;
	for (i = 0; i < NR; i++) {
	  mr = MR_New(mesh);
	  
	  status = fscanf(fp,"%d",&nrv);
	  if (status == EOF)
	    MSTK_Report("MESH_InitFromFile",
			"Premature end of file while reading region data"
			,FATAL);
	  else if (status == 0)
	    MSTK_Report("MESH_InitFromFile",
			"Error in reading region data",FATAL);

	  if (rverts) {
	    if (nrv > max_nrv) {
	      max_nrv = nrv;
	      rverts = MSTK_realloc(rverts,max_nrv*sizeof(MVertex_ptr));
	    }
	  }
	  else {
	    max_nrv = nrv;
	    rverts = MSTK_malloc(nrv*sizeof(MVertex_ptr));
	  }
	  
	  for (j = 0; j < nrv; j++) {
	    status = fscanf(fp,"%d",&vid1);
	    if (status == EOF)
	      MSTK_Report("MESH_InitFromFile",
		      "Premature end of file while reading region data",
			  FATAL);
	    else if (status == 0)
	      MSTK_Report("MESH_InitFromFile",
			  "Error in reading region data",FATAL);

	    rverts[j] = List_Entry(mesh->mvertex,vid1-1);

	  }
	  MR_Set_Vertices(mr,nrv,rverts,0,NULL);  
	  
	  status = fscanf(fp,"%d %d %s %d %d %d",&gdim,&gid,
			  ptype_str,&id,&masterparid,&globalid);
	  if (status == EOF)
	    MSTK_Report("MESH_InitFromFile",
			"Premature end of file while reading regions",FATAL);
	  else if (status == 0)
	    MSTK_Report("MESH_InitFromFile",
			"Error in reading region data",FATAL);
	  
	  MR_Set_GEntDim(mr,gdim);
	  MR_Set_ID(mr,id);
	  MR_Set_PType(mr,Str2PType(ptype_str));
	  MR_Set_MasterParID(mr,masterparid);
	  MR_Set_GlobalID(mr,globalid);
	  MR_Set_GEntID(mr,gid);
	  
	}
	if (rverts)
	  MSTK_free(rverts);

	processed_regions = 1;
      }
      else {
	MSTK_Report("MESH_InitFromFile",
		    "Expected region description in terms of faces",ERROR);
        fclose(fp);
	return 0;
      }
    }
    else if (strncmp(rltype_str,"face",4) == 0) {
      if (mesh->reptype != R1 && mesh->reptype != R2) {
	rfaces = NULL;
	rfdirs = NULL;
	for (i = 0; i < NR; i++) {
	  mr = MR_New(mesh);
	  
	  status = fscanf(fp,"%d",&nrf);
	  if (status == EOF)
	    MSTK_Report("MESH_InitFromFile",
			"Premature end of file while reading region data",
			FATAL);
	  else if (status == 0)
	    MSTK_Report("MESH_InitFromFile",
			"Error in reading region data",FATAL);

	  if (rfaces) {
	    if (nrf > max_nrf) {
	      max_nrf = nrf;
	      rfaces = MSTK_realloc(rfaces,max_nrf*sizeof(MFace_ptr));
	      rfdirs = MSTK_realloc(rfdirs,max_nrf*sizeof(int));
	    }
	  }
	  else {
	    max_nrf = nrf;
	    rfaces = MSTK_malloc(nrf*sizeof(MFace_ptr));
	    rfdirs = MSTK_malloc(nrf*sizeof(int));
	  }
	
	  for (j = 0; j < nrf; j++) {
	    status = fscanf(fp,"%d",&fid);
	    if (status == EOF)
	      MSTK_Report("MESH_InitFromFile",
			  "Premature end of file while reading region data",
			  FATAL);
	    else if (status == 0)
	      MSTK_Report("MESH_InitFromFile",
			  "Error in reading region data",FATAL);

	    rfdirs[j] = fid > 0 ? 1 : 0;
	    rfaces[j] = List_Entry(mesh->mface,abs(fid)-1);
	  }
	  
	  MR_Set_Faces(mr,nrf,rfaces,rfdirs);
	  
	  status = fscanf(fp,"%d %d %s %d %d %d",&gdim,&gid,
		  ptype_str,&id,&masterparid,&globalid);
	  if (status == EOF)
	    MSTK_Report("MESH_InitFromFile",
			"Premature end of file while reading faces",FATAL);
	  else if (status == 0)
	    MSTK_Report("MESH_InitFromFile",
			"Error in reading region data",FATAL);
	  
	  MR_Set_GEntDim(mr,gdim);
	  MR_Set_ID(mr,id);
	  MR_Set_PType(mr,Str2PType(ptype_str));
	  MR_Set_MasterParID(mr,masterparid);
	  MR_Set_GlobalID(mr,globalid);
	 
	  MR_Set_GEntID(mr,gid);
	}
	if (rfaces) {
	  MSTK_free(rfaces);
	  MSTK_free(rfdirs);
	}

	processed_regions = 1;
      }
    }
    else {
      MSTK_Report("MESH_InitFromFile",
		  "Expected region description in terms of vertices",ERROR);
      fclose(fp);
      return 0;
    }
  }

  if (processed_regions) {
    status = fscanf(fp,"%s",temp_str);
    if (status == EOF) {
      if (mesh->reptype == R2)
	MSTK_Report("MESH_InitFromFile",
		    "Premature end of file after reading regions",FATAL);
      else {
        fclose(fp);
	return 1;
      }
    }
  }


  /* ADJACENT REGION DATA */

  if (strncmp(temp_str,"adjregions",10) == 0) {
    if (mesh->reptype == R2) {
      for (i = 0; i < NR; i++) {
	mr = List_Entry(mesh->mregion,i);

	status = fscanf(fp,"%d",&nar);
	if (status == EOF)
	  MSTK_Report("MESH_InitFromFile",
		      "Premature end of file while reading adjacent regions",
		      FATAL);
	else if (status == 0)
	  MSTK_Report("MESH_InitFromFile",
		      "Error in reading adjacent region data",FATAL);

	for (j = 0; j < nar; j++) {
	  status = fscanf(fp,"%d",&adjrid);
	  if (status == EOF)
	    MSTK_Report("MESH_InitFromFile",
			"Premature end of file while reading ajdacent regions",
			FATAL);
	  else if (status == 0)
	    MSTK_Report("MESH_InitFromFile",
			"Error in reading adjacent region data",FATAL);
	  
	  if (adjrid == 0)
	    continue;
	  adjr = List_Entry(mesh->mregion,adjrid-1);
	  MR_Add_AdjRegion(mr,j,adjr);
	}
      }

      processed_adjr = 1;
    }
    else {
      for (i = 0; i < NR; i++) {
	status = fscanf(fp,"%d",&nar);
	if (status == EOF)
	  MSTK_Report("MESH_InitFromFile",
		      "Premature end of file while reading adjacent regions",
		      FATAL);
	else if (status == 0)
	  MSTK_Report("MESH_InitFromFile",
		      "Error in reading adjacent region data",FATAL);
	
	for (j = 0; j < nar; j++) {
	  fscanf(fp,"%d",&rid);
	  if (status == EOF)
	    MSTK_Report("MESH_InitFromFile",
			"Premature end of file while reading adjacent regions",
			FATAL);
	  else if (status == 0)
	    MSTK_Report("MESH_InitFromFile",
			"Error in reading adjacent region data",FATAL);
	}
      }

      processed_adjr = 1;
    }
  }

  if (processed_adjr) {
    status = fscanf(fp,"%s",temp_str);
    if (status == EOF) {
      if (mesh->reptype == R2)
	MSTK_Report("MESH_InitFromFile",
		    "Premature end of file after reading regions",FATAL);
      else {
        fclose(fp);
	return 1;
      }
    }
  }

  /* ATTRIBUTE DATA */

  if (strncmp(temp_str,"attributes",10) == 0) {

    done = 0;
    while (!done) {
      status = fscanf(fp,"%s",attname);
      if (status == EOF) {
	done = 1;
	continue;
      }
      else if (status == 0)
	MSTK_Report("MESH_InitFromFile",
		    "Error in reading attribute data",FATAL);

      status = fscanf(fp,"%s",atttype_str);
      if (status == EOF)
	MSTK_Report("MESH_InitFromFile",
		    "Premature end of file while reading attributes",FATAL);
      else if (status == 0)
	MSTK_Report("MESH_InitFromFile",
		    "Error in reading attribute data",FATAL);

      if (strncmp(atttype_str,"INT",3) == 0)
	atttype = INT;
      else if (strncmp(atttype_str,"DOUBLE",6) == 0)
	atttype = DOUBLE;
      else if (strncmp(atttype_str,"POINTER",7) == 0)
	MSTK_Report("MESH_InitFromFile",
		    "Cannot specify POINTER attributes in file",FATAL);
      else if (strncmp(atttype_str,"VECTOR",6) == 0)
	atttype = VECTOR;
      else if (strncmp(atttype_str,"TENSOR",6) == 0)
	atttype = TENSOR;
      else {
	sprintf(temp_str,"%-s not a recognized attribute type",atttype_str);
	MSTK_Report("MESH_InitFromFile",temp_str,FATAL);
      }

      status = fscanf(fp,"%d",&ncomp);
      if (status == EOF)
	MSTK_Report("MESH_InitFromFile",
		    "Premature end of file while reading attributes",FATAL);
      else if (status == 0)
	MSTK_Report("MESH_InitFromFile",
		    "Error in reading attribute data",FATAL);

      if ((atttype == INT || atttype == DOUBLE) && (ncomp != 1)) 
	MSTK_Report("MESH_InitFromFile","Number of components should be 1 for attributes of type INT or DOUBLE",WARN);
      else if ((atttype == VECTOR || atttype == TENSOR) && (ncomp == 0))
	MSTK_Report("MESH_InitFromFile","Number of components should be non-zero for attributes of type VECTOR or TENSOR",FATAL);



      status = fscanf(fp,"%s",attent_str);
      if (status == EOF)
	MSTK_Report("MESH_InitFromFile",
		    "Premature end of file while reading attributes",FATAL);
      else if (status == 0)
	MSTK_Report("MESH_InitFromFile",
		    "Error in reading attribute data",FATAL);
      if (strncmp(attent_str,"MVERTEX",7) == 0) 
	attent = MVERTEX;
      else if (strncmp(attent_str,"MEDGE",5) == 0)
	attent = MEDGE;
      else if (strncmp(attent_str,"MFACE",5) == 0)
	attent = MFACE;
      else if (strncmp(attent_str,"MREGION",7) == 0)
	attent = MREGION;
      else if (strncmp(attent_str,"MALLTYPE",8) == 0)
	attent = MALLTYPE;
      else {
	sprintf(temp_str,"%s not a recognized entity type for attributes",
		attent_str);
	MSTK_Report("MESH_InitFromFile",temp_str,FATAL);
      }

      status = fscanf(fp,"%d",&nent);
      if (status == EOF)
	MSTK_Report("MESH_InitFromFile",
		    "Premature end of file while reading attributes",FATAL);
      else if (status == 0)
	MSTK_Report("MESH_InitFromFile",
		    "Error in reading attribute data",FATAL);

      if (nent < 1) 
	MSTK_Report("MESH_InitFromFile",
		    "Attribute applied on no entities?",FATAL);

      if (atttype == INT || atttype == DOUBLE || atttype == POINTER)
	attrib = MAttrib_New(mesh,attname,atttype,attent);
      else
	attrib = MAttrib_New(mesh,attname,atttype,attent,ncomp);
      
      for (i = 0; i < nent; i++) {
	fscanf(fp,"%d %d",&dim,&id);
	if (status == EOF)
	  MSTK_Report("MESH_InitFromFile",
		      "Premature end of file while reading attributes",FATAL);
	else if (status == 0)
	  MSTK_Report("MESH_InitFromFile",
		      "Error in reading attribute data",FATAL);
	
	if (attent != dim && attent != MALLTYPE) {
	  MSTK_Report("MESH_InitFromFile",
		      "Attribute not applicable to this type of entity",WARN);
	  if (atttype == INT)
	    fscanf(fp,"%d",&ival);
	  else
	    for (j = 0; j < ncomp; j++)
	      fscanf(fp,"%lf",&rval);
	}
	else {

	  if (atttype == INT)
	    fscanf(fp,"%d",&ival);
	  else if (atttype == DOUBLE)
	    fscanf(fp,"%lf",&rval);
	  else if (atttype == VECTOR || atttype == TENSOR) {
	    ival = ncomp;
	    rval_arr = (double *) MSTK_malloc(ncomp*sizeof(double));
	    for (j = 0; j < ncomp; j++) {
	      fscanf(fp,"%lf",&(rval_arr[j]));
	    }
	  }

	  switch (dim) {
	  case MVERTEX:
	    ent = MESH_VertexFromID(mesh,id);
	    break;
	  case MEDGE:
	    ent = MESH_EdgeFromID(mesh,id);
	    break;
	  case MFACE:
	    ent = MESH_FaceFromID(mesh,id);
	    break;
	  case MREGION:
	    ent = MESH_RegionFromID(mesh,id);
	    break;
	  default:
	    ent = NULL;
	    MSTK_Report("MESH_InitFromFile","Invalid entity type",FATAL);
	  }
	  
	  MEnt_Set_AttVal(ent,attrib,ival,rval,(void *)rval_arr);
	  
	}

      } /* for (i = 0; i < nent; i++) */

    } /* while (!done) */

  }

  fclose(fp);

  return 1;
}

int PMESH_Ghost_WriteToFile(Mesh_ptr mesh, const char *filename, RepType rtype) {
  FILE *fp;
  char mesg[80], attname[256];
  int i, j, k, idx;
  int gdim, gid;
  int mvid, mvid0, mvid1, mvid2, mrid2, meid, mfid, mrid;
  int nav, nar, nfe, nfv, nrf, nrv, dir=0;
  int nv, ne, nf, nr;
  int natt, ncomp, ival, nent;
  double xyz[3], rval, rdummy, *rval_arr;
  void *pval, *pdummy;
  MVertex_ptr mv, mv0, mv1, mv2;
  MEdge_ptr me;
  MFace_ptr mf;
  MRegion_ptr mr, mr2;
  GEntity_ptr gent;
  List_ptr adjverts, mfedges, mfverts, mrfaces, mrverts, adjregs;
  RepType reptype;
  MAttrib_ptr attrib, vidatt, eidatt, fidatt, ridatt;
  MType attentdim;
  MAttType atttype;
  /*for mpi */
  PType ptype;
  int id, masterparid, globalid;
  /*end for mpi */
  

  if (!(fp = fopen(filename,"w"))) {
    sprintf(mesg,"Cannot open file %-s for writing",filename);
    MSTK_Report("MESH_WriteToFile",mesg,ERROR);
    return 0;
  }

  if (rtype != UNKNOWN_REP) {
    reptype = rtype;
  }
  else {
    reptype = MESH_RepType(mesh);
  }

  nv = MESH_Num_Vertices(mesh);
  ne = MESH_Num_Edges(mesh);
  nf = MESH_Num_Faces(mesh);
  nr = MESH_Num_Regions(mesh);

  fprintf(fp,"MSTK %-2.1lf\n",MSTK_VER);
  fprintf(fp,"%s %d %d %d %d\n",
	  MESH_rtype_str[reptype], 
	  nv+MESH_Num_GhostVertices(mesh),
	  (reptype >= R1 && reptype <= R4)?0:ne+MESH_Num_GhostEdges(mesh), 
	  (reptype >= R1 && reptype <= R2 && nr)?0:nf+MESH_Num_GhostFaces(mesh), 
	  nr+MESH_Num_GhostRegions(mesh));
  vidatt = MAttrib_New(mesh,"vidatt",INT,MVERTEX);
  eidatt = MAttrib_New(mesh,"eidatt",INT,MEDGE);
  fidatt = MAttrib_New(mesh,"fidatt",INT,MFACE);
  ridatt = MAttrib_New(mesh,"ridatt",INT,MREGION);

  idx = 0; i = 0;
  while ((mv = MESH_Next_Vertex(mesh,&idx)))
    MEnt_Set_AttVal(mv,vidatt,++i,0.0,NULL);
  /* for ghost */
  idx = 0; i = 0;
  while ((mv = MESH_Next_GhostVertex(mesh,&idx)))
    MEnt_Set_AttVal(mv,vidatt,++i,0.0,NULL);

  idx = 0; i = 0;
  while ((me = MESH_Next_Edge(mesh,&idx)))
    MEnt_Set_AttVal(me,eidatt,++i,0.0,NULL);
  /* for ghost */
  idx = 0; i = 0;
  while ((me = MESH_Next_GhostEdge(mesh,&idx)))
    MEnt_Set_AttVal(me,eidatt,++i,0.0,NULL);

  idx = 0; i = 0;
  while ((mf = MESH_Next_Face(mesh,&idx)))
    MEnt_Set_AttVal(mf,fidatt,++i,0.0,NULL);
  /* for ghost */
  idx = 0; i = 0;
  while ((mf = MESH_Next_GhostFace(mesh,&idx)))
    MEnt_Set_AttVal(mf,fidatt,++i,0.0,NULL);

  idx = 0; i = 0;
  while ((mr = MESH_Next_Region(mesh,&idx)))
    MEnt_Set_AttVal(mr,ridatt,++i,0.0,NULL);
  /* for ghost */
  idx = 0; i = 0;
  while ((mr = MESH_Next_GhostRegion(mesh,&idx)))
    MEnt_Set_AttVal(mr,ridatt,++i,0.0,NULL);
  
  fprintf(fp,"vertices\n");
  idx = 0;
  while ((mv = MESH_Next_Vertex(mesh,&idx))) {

    MV_Coords(mv,xyz);

    gdim = MV_GEntDim(mv);
    gid = MV_GEntID(mv);
    /* for mpi */
    
    ptype = MV_PType(mv);
    id = MV_ID(mv);
    masterparid = MV_MasterParID(mv);
    globalid = MV_GlobalID(mv);
    fprintf(fp,"%24.16lf %24.16lf %24.16lf   %d %d\t%12s %6d %3d %6d\n",
	    xyz[0],xyz[1],xyz[2],gdim,gid,MESH_ptype_str[ptype],id,masterparid,globalid);
    /* end for mpi */
    /*    fprintf(fp,"%24.16lf %24.16lf %24.16lf   %d %d\n",
	  xyz[0],xyz[1],xyz[2],gdim,gid); */
    
  }
  /* for ghost */
  idx = 0;
  while ((mv = MESH_Next_GhostVertex(mesh,&idx))) {

    MV_Coords(mv,xyz);

    gdim = MV_GEntDim(mv);
    gid = MV_GEntID(mv);
    /* for mpi */
    
    ptype = MV_PType(mv);
    id = MV_ID(mv);
    masterparid = MV_MasterParID(mv);
    globalid = MV_GlobalID(mv);
    fprintf(fp,"%24.16lf %24.16lf %24.16lf   %d %d\t%12s %6d %3d %6d\n",
	    xyz[0],xyz[1],xyz[2],gdim,gid,MESH_ptype_str[ptype],id,masterparid,globalid);
    /* end for mpi */
    /*    fprintf(fp,"%24.16lf %24.16lf %24.16lf   %d %d\n",
	  xyz[0],xyz[1],xyz[2],gdim,gid); */
    
  }

  if (reptype == R2 || reptype == R4) {
    fprintf(fp,"adjvertices\n");

    idx = 0;
    while ((mv = MESH_Next_Vertex(mesh,&idx))) {

      nav = MV_Num_AdjVertices(mv);
      fprintf(fp,"%d ",nav);
      
      adjverts = MV_AdjVertices(mv);
      for (j = 0; j < nav; j++) {
	mv2 = List_Entry(adjverts,j);
	MEnt_Get_AttVal(mv2,vidatt,&mvid2,&rval,&pval);
	fprintf(fp,"%d ",mvid2);
      }
      fprintf(fp,"\n");
      List_Delete(adjverts);
    }
  }



  if (reptype <= F4 && ne) {
    fprintf(fp,"edges\n");

    idx = 0;
    while ((me = MESH_Next_Edge(mesh,&idx))) {

      mv0 = ME_Vertex(me,0);
      MEnt_Get_AttVal(mv0,vidatt,&mvid0,&rval,&pval);
      mv1 = ME_Vertex(me,1);
      MEnt_Get_AttVal(mv1,vidatt,&mvid1,&rval,&pval);

      gdim = ME_GEntDim(me);
      gid = ME_GEntID(me);
      /* for mpi */
      ptype = ME_PType(me);
      id = ME_ID(me);
      masterparid = ME_MasterParID(me);
      globalid = ME_GlobalID(me);
      fprintf(fp,"%d %d \t%d %d\t%12s %6d %3d %6d\n"
	      ,mvid0,mvid1,gdim,gid,MESH_ptype_str[ptype],id,masterparid,globalid);
      /* end for mpi */
      /* fprintf(fp,"%d %d \t%d %d\n",mvid0,mvid1,gdim,gid); */
    }
  }



  if (reptype <= F4) {

    /* For full representations, always write out faces in terms of edges */

    fprintf(fp,"faces edge\n");
    
    idx = 0;
    while ((mf = MESH_Next_Face(mesh,&idx))) {
      
      nfe = MF_Num_Edges(mf);
      fprintf(fp,"%d ",nfe);
      
      mfedges = MF_Edges(mf,1,0);
      for (j = 0; j < nfe; j++) {
	me = List_Entry(mfedges,j);
	dir = MF_EdgeDir_i(mf,j);
	MEnt_Get_AttVal(me,eidatt,&meid,&rval,&pval);
	if (dir != 1) meid = -meid;
	fprintf(fp,"%d ",meid);
      }
      List_Delete(mfedges);
      
      gdim = MF_GEntDim(mf);
      /*
	gent = MF_GEntity(mf);
	gid = gent ? -99 : 0;
      */
      gid = MF_GEntID(mf);

      /* for mpi */
      ptype = MF_PType(mf);
      id = MF_ID(mf);
      masterparid = MF_MasterParID(mf);
      globalid = MF_GlobalID(mf);
      fprintf(fp,"\t%d %d\t%12s %6d %3d %6d\n",
	      gdim,gid,MESH_ptype_str[ptype],id,masterparid,globalid);
      /* end for mpi */
      /* fprintf(fp,"\t%d %d\n",gdim,gid); */
    }
    idx = 0;
    while ((mf = MESH_Next_GhostFace(mesh,&idx))) {
      
      nfe = MF_Num_Edges(mf);
      fprintf(fp,"%d ",nfe);
      
      mfedges = MF_Edges(mf,1,0);
      for (j = 0; j < nfe; j++) {
	me = List_Entry(mfedges,j);
	dir = MF_EdgeDir_i(mf,j);
	MEnt_Get_AttVal(me,eidatt,&meid,&rval,&pval);
	if (dir != 1) meid = -meid;
	fprintf(fp,"%d ",meid);
      }
      List_Delete(mfedges);
      
      gdim = MF_GEntDim(mf);
      /*
	gent = MF_GEntity(mf);
	gid = gent ? -99 : 0;
      */
      gid = MF_GEntID(mf);

      /* for mpi */
      ptype = MF_PType(mf);
      id = MF_ID(mf);
      masterparid = MF_MasterParID(mf);
      globalid = MF_GlobalID(mf);
      fprintf(fp,"\t%d %d\t%12s %6d %3d %6d\n",
	      gdim,gid,MESH_ptype_str[ptype],id,masterparid,globalid);
      /* end for mpi */
      /* fprintf(fp,"\t%d %d\n",gdim,gid); */
    }
  }
  else {

    /* For reduced representations, R3 and R4 always write out faces
       in terms of vertices. For reduced representations, R1 and R2
       write out faces in terms of vertices only when there are no
       regions (i.e. faces are the highest level mesh entities) */

    if ((reptype > R2) || (nr == 0)) {

      fprintf(fp,"faces vertex\n");

      idx = 0;
      while ((mf = MESH_Next_Face(mesh,&idx))) {
	
	nfv = MF_Num_Edges(mf);
	fprintf(fp,"%d ",nfv);
	
	mfverts = MF_Vertices(mf,1,0);
	for (j = 0; j < nfv; j++) {
	  mv = List_Entry(mfverts,j);
	  MEnt_Get_AttVal(mv,vidatt,&mvid,&rval,&pval);
	  fprintf(fp,"%d ",mvid);
	}
	List_Delete(mfverts);

	gdim = MF_GEntDim(mf);
	gid = MF_GEntID(mf);
	/* for mpi */
	ptype = MF_PType(mf);
	id = MF_ID(mf);
	masterparid = MF_MasterParID(mf);
	globalid = MF_GlobalID(mf);
	fprintf(fp,"\t%d %d\t%12s %6d %3d %6d\n",
		gdim,gid,MESH_ptype_str[ptype],id,masterparid,globalid);
	/* end for mpi */
	/* fprintf(fp,"\t%d %d\n",gdim,gid);*/
      }
    idx = 0;
    while ((mf = MESH_Next_GhostFace(mesh,&idx))) {
      
      nfv = MF_Num_Edges(mf);
      fprintf(fp,"%d ",nfv);
      
      mfverts = MF_Vertices(mf,1,0);
      for (j = 0; j < nfv; j++) {
	mv = List_Entry(mfverts,j);
	MEnt_Get_AttVal(mv,vidatt,&mvid,&rval,&pval);
	fprintf(fp,"%d ",mvid);
      }
      List_Delete(mfverts);
      
      gdim = MF_GEntDim(mf);
      gid = MF_GEntID(mf);
      /* for mpi */
      ptype = MF_PType(mf);
      id = MF_ID(mf);
      masterparid = MF_MasterParID(mf);
      globalid = MF_GlobalID(mf);
      fprintf(fp,"\t%d %d\t%12s %6d %3d %6d\n",
	      gdim,gid,MESH_ptype_str[ptype],id,masterparid,globalid);
      /* end for mpi */
      /* fprintf(fp,"\t%d %d\n",gdim,gid);*/
    }
    }
  }
  
  if (nr) {
    if (reptype <= F4 || reptype >= R2) {
      fprintf(fp,"regions face\n");

      idx = 0;
      while ((mr = MESH_Next_Region(mesh,&idx))) {

	nrf = MR_Num_Faces(mr);
	fprintf(fp,"%d ",nrf);

	mrfaces = MR_Faces(mr);
	for (j = 0; j < nrf; j++) {
	  mf = List_Entry(mrfaces,j);
	  dir = MR_FaceDir_i(mr,j);
	  MEnt_Get_AttVal(mf,fidatt,&mfid,&rval,&pval);
	  if (dir != 1) mfid = -mfid;
	  fprintf(fp,"%d ",mfid);
	}
	List_Delete(mrfaces);
	
	gdim = MF_GEntDim(mr);
	gid = MR_GEntID(mr);
	/* for mpi */
	ptype = MR_PType(mr);
	id = MR_ID(mr);
	masterparid = MR_MasterParID(mr);
	globalid = MR_GlobalID(mr);
	fprintf(fp,"\t%d %d\t%12s %6d %3d %6d\n",
		gdim,gid,MESH_ptype_str[ptype],id,masterparid,globalid);
	/* end for mpi */
	/* fprintf(fp,"\t%d %d\n",gdim,gid); */
      }
    }
    else {
      fprintf(fp,"regions vertex\n");

      idx = 0;
      while ((mr = MESH_Next_Region(mesh,&idx))) {

	nrv = MR_Num_Vertices(mr);
	fprintf(fp,"%d ",nrv);

	mrverts = MR_Vertices(mr);
	for (j = 0; j < nrv; j++) {
	  mv = List_Entry(mrverts,j);
	  MEnt_Get_AttVal(mv,vidatt,&mvid,&rval,&pval);
	  fprintf(fp,"%d ",mvid);
	}
	List_Delete(mrverts);
	
	gdim = MR_GEntDim(mr);
	gid = MR_GEntID(mr);
	/* for mpi */
	ptype = MR_PType(mr);
	id = MR_ID(mr);
	masterparid = MR_MasterParID(mr);
	globalid = MR_GlobalID(mr);
	fprintf(fp,"\t%d %d\t%12s %6d %3d %6d\n",
		gdim,gid,MESH_ptype_str[ptype],id,masterparid,globalid);
	/* end for mpi */

	/*fprintf(fp,"\t%d %d\n",gdim,gid);*/
      }
    }

    if (reptype == R2 || reptype == R4) {
      fprintf(fp,"adjregions\n");
      
      idx = 0;
      while ((mr = MESH_Next_Region(mesh,&idx))) {

	nar = MR_Num_Faces(mr);
	fprintf(fp,"%d ",nar);

	adjregs = MR_AdjRegions(mr);

	for (j = 0; j < nar; j++) {
	  mr2 = List_Entry(adjregs,j);
	  if ((long) mr2 == -1) 
	    fprintf(fp,"%d ",0);
	  else {
	    MEnt_Get_AttVal(mr2,ridatt,&mrid2,&rval,&pval);
	    fprintf(fp,"%d ",mrid2);
	  }
	}
	fprintf(fp,"\n");
	List_Delete(adjregs);
      }
    }
  }


  /* Write out attributes if there are more than the 4 that we created 
    in this routine */


  if ((natt = MESH_Num_Attribs(mesh)) > 4) {

    fprintf(fp,"attributes\n");

    for (i = 0; i < natt; i++) {
      
      attrib = MESH_Attrib(mesh,i);

      /* Don't write out attribs we created for the internal use of 
	 this routine */
      if (attrib == vidatt || attrib == eidatt || attrib == fidatt || 
	  attrib == ridatt) continue;
      
      MAttrib_Get_Name(attrib,attname);

      atttype = MAttrib_Get_Type(attrib);
      if (atttype == POINTER) continue;  /* cannot write it out */

      ncomp = MAttrib_Get_NumComps(attrib);

      attentdim = MAttrib_Get_EntDim(attrib);


      /* First count how many entities actually have the attribute assigned */

      nent = 0;
      switch(attentdim) {
      case MVERTEX:
	idx = 0;
	while ((mv = MESH_Next_Vertex(mesh,&idx)))
	  if (MEnt_Get_AttVal(mv,attrib,&ival,&rval,&pval)) nent++;
	break;
      case MEDGE:
	idx = 0;
	while ((me = MESH_Next_Edge(mesh,&idx)))
	  if (MEnt_Get_AttVal(me,attrib,&ival,&rval,&pval)) nent++;
	break;
      case MFACE:
	idx = 0;
	while ((mf = MESH_Next_Face(mesh,&idx)))
	  if (MEnt_Get_AttVal(mf,attrib,&ival,&rval,&pval)) nent++;	    
	break;
      case MREGION: 
	idx = 0;
	while ((mr = MESH_Next_Region(mesh,&idx)))
	  if (MEnt_Get_AttVal(mr,attrib,&ival,&rval,&pval)) nent++;
	break;
      case MALLTYPE:
	idx = 0;
	while ((mv = MESH_Next_Vertex(mesh,&idx)))
	  if (MEnt_Get_AttVal(mv,attrib,&ival,&rval,&pval)) nent++;
	idx = 0;
	while ((me = MESH_Next_Edge(mesh,&idx)))
	  if (MEnt_Get_AttVal(me,attrib,&ival,&rval,&pval)) nent++;
	idx = 0;
	while ((mf = MESH_Next_Face(mesh,&idx)))
	  if (MEnt_Get_AttVal(mf,attrib,&ival,&rval,&pval)) nent++;	    
	idx = 0;
	while ((mr = MESH_Next_Region(mesh,&idx)))
	  if (MEnt_Get_AttVal(mr,attrib,&ival,&rval,&pval)) nent++;
	break;	
      default:
	break;
      } /* switch (attentdim) */


      /* No point in writing out attribute if no entity uses it! Or is there? */

      if (!nent) continue;



      fprintf(fp,"%-s\n",attname);

      switch(atttype) {
      case INT:
	fprintf(fp,"INT\n");
	break;
      case DOUBLE:
	fprintf(fp,"DOUBLE\n");
	break;
      case VECTOR:
	fprintf(fp,"VECTOR\n");
	break;
      case TENSOR:
	fprintf(fp,"TENSOR\n");
	break;
      default:
	MSTK_Report("MESH_WriteToFile",
		    "Unrecognizable or unprintable attribute type\n",WARN);
	continue;	
      }

      fprintf(fp,"%-d\n",ncomp);

      switch(attentdim) {
      case MVERTEX:
	fprintf(fp,"MVERTEX\n");
	break;
      case MEDGE:
	fprintf(fp,"MEDGE\n");
	break;
      case MFACE:
	fprintf(fp,"MFACE\n");
	break;
      case MREGION:
	fprintf(fp,"MREGION\n");
	break;
      case MALLTYPE:
	fprintf(fp,"MALLTYPE\n");
	break;
      default:
	MSTK_Report("Mesh_WriteToFile","Unrecognized entity type",WARN);
	break;
      }

      fprintf(fp,"%-d\n",nent);


      switch(attentdim) {
      case MVERTEX:
	idx = 0;
	while ((mv = MESH_Next_Vertex(mesh,&idx))) {
	  if (MEnt_Get_AttVal(mv,attrib,&ival,&rval,&pval)) {
	    MEnt_Get_AttVal(mv,vidatt,&mvid,&rdummy,&pdummy);
	    fprintf(fp,"0 %-d ",mvid);
	    switch (atttype) {
	    case INT:
	      fprintf(fp," %-d",ival);
	      break;
	    case DOUBLE: 
	      fprintf(fp," %-lf ",rval);
	      break;
	    case VECTOR: case TENSOR:
	      rval_arr = (double *) pval;
	      for (k = 0; k < ncomp; k++)
		fprintf(fp," %-lf ",rval_arr[k]);
	      break;
	    default:
	      break;
	    }
	    fprintf(fp,"\n");
	  }
	}
	break;
      case MEDGE:
	idx = 0;
	while ((me = MESH_Next_Edge(mesh,&idx))) {
	  if (MEnt_Get_AttVal(me,attrib,&ival,&rval,&pval)) {
	    MEnt_Get_AttVal(me,eidatt,&meid,&rdummy,&pdummy);
	    fprintf(fp,"1 %-d ",meid);
	    switch (atttype) {
	    case INT:
	      fprintf(fp," %-d",ival);
	      break;
	    case DOUBLE: 
	      fprintf(fp," %-lf ",rval);
	      break;
	    case VECTOR: case TENSOR:
	      rval_arr = (double *) pval;
	      for (k = 0; k < ncomp; k++)
		fprintf(fp," %-lf ",rval_arr[k]);
	      break;
	    default:
	      break;
	    }
	    fprintf(fp,"\n");
	  }
	}
	break;
      case MFACE:
	idx = 0;
	while ((mf = MESH_Next_Face(mesh,&idx))) {
	  if (MEnt_Get_AttVal(mf,attrib,&ival,&rval,&pval)) {
	    MEnt_Get_AttVal(mf,fidatt,&mfid,&rdummy,&pdummy);
	    fprintf(fp,"2 %-d ",mfid);
	    switch (atttype) {
	    case INT:
	      fprintf(fp," %-d",ival);
	      break;
	    case DOUBLE: 
	      fprintf(fp," %-lf ",rval);
	      break;
	    case VECTOR: case TENSOR:
	      rval_arr = (double *) pval;
	      for (k = 0; k < ncomp; k++)
		fprintf(fp," %-lf ",rval_arr[k]);
	      break;
	    default:
	      break;
	    }
	    fprintf(fp,"\n");
	  }
	}
	break;

      case MREGION: 
	idx = 0;
	while ((mr = MESH_Next_Region(mesh,&idx))) {
	  if (MEnt_Get_AttVal(mr,attrib,&ival,&rval,&pval)) {
	    MEnt_Get_AttVal(mr,ridatt,&mrid,&rdummy,&pdummy);
	    fprintf(fp,"3 %-d ",mrid);
	    switch (atttype) {
	    case INT:
	      fprintf(fp," %-d",ival);
	      break;
	    case DOUBLE: 
	      fprintf(fp," %-lf ",rval);
	      break;
	    case VECTOR: case TENSOR:
	      rval_arr = (double *) pval;
	      for (k = 0; k < ncomp; k++)
		fprintf(fp," %-lf ",rval_arr[k]);
	      break;
	    default:
	      break;
	    }
	    fprintf(fp,"\n");
	  }
	}
	break;

      case MALLTYPE:
	idx = 0;
	while ((mv = MESH_Next_Vertex(mesh,&idx))) {
	  if (MEnt_Get_AttVal(mv,attrib,&ival,&rval,&pval)) {
	    MEnt_Get_AttVal(mv,vidatt,&mvid,&rdummy,&pdummy);
	    fprintf(fp,"0 %-d ",mvid);
	    switch (atttype) {
	    case INT:
	      fprintf(fp," %-d",ival);
	      break;
	    case DOUBLE: 
	      fprintf(fp," %-lf ",rval);
	      break;
	    case VECTOR: case TENSOR:
	      rval_arr = (double *) pval;
	      for (k = 0; k < ncomp; k++)
		fprintf(fp," %-lf ",rval_arr[k]);
	      break;
	    default:
	      break;
	    }
	    fprintf(fp,"\n");
	  }
	}
	idx = 0;
	while ((me = MESH_Next_Edge(mesh,&idx))) {
	  if (MEnt_Get_AttVal(me,attrib,&ival,&rval,&pval)) {
	    MEnt_Get_AttVal(me,eidatt,&meid,&rdummy,&pdummy);
	    fprintf(fp,"1 %-d ",meid);
	    switch (atttype) {
	    case INT:
	      fprintf(fp," %-d",ival);
	      break;
	    case DOUBLE: 
	      fprintf(fp," %-lf ",rval);
	      break;
	    case VECTOR: case TENSOR:
	      rval_arr = (double *) pval;
	      for (k = 0; k < ncomp; k++)
		fprintf(fp," %-lf ",rval_arr[k]);
	      break;
	    default:
	      break;
	    }
	    fprintf(fp,"\n");
	  }
	}
	idx = 0;
	while ((mf = MESH_Next_Face(mesh,&idx))) {
	  if (MEnt_Get_AttVal(mf,attrib,&ival,&rval,&pval)) {
	    MEnt_Get_AttVal(mf,fidatt,&mfid,&rdummy,&pdummy);
	    fprintf(fp,"2 %-d ",mfid);
	    switch (atttype) {
	    case INT:
	      fprintf(fp," %-d",ival);
	      break;
	    case DOUBLE: 
	      fprintf(fp," %-lf ",rval);
	      break;
	    case VECTOR: case TENSOR:
	      rval_arr = (double *) pval;
	      for (k = 0; k < ncomp; k++)
		fprintf(fp," %-lf ",rval_arr[k]);
	      break;
	    default:
	      break;
	    }
	    fprintf(fp,"\n");
	  }
	}
	idx = 0;
	while ((mr = MESH_Next_Region(mesh,&idx))) {
	  if (MEnt_Get_AttVal(mr,attrib,&ival,&rval,&pval)) {
	    MEnt_Get_AttVal(mr,ridatt,&mrid,&rdummy,&pdummy);
	    fprintf(fp,"3 %-d ",mrid);
	    switch (atttype) {
	    case INT:
	      fprintf(fp," %-d",ival);
	      break;
	    case DOUBLE: 
	      fprintf(fp," %-lf ",rval);
	      break;
	    case VECTOR: case TENSOR:
	      rval_arr = (double *) pval;
	      for (k = 0; k < ncomp; k++)
		fprintf(fp," %-lf ",rval_arr[k]);
	      break;
	    default:
	      break;
	    }
	    fprintf(fp,"\n");
	  }
	}
	break;	
      default:
	break;
      } /* switch (attentdim) */

    } /* for (i = 0; i < natt) */
    
  } /* if (Mesh_Num_Attribs(mesh)) */
  

  idx = 0; i = 0;
  while ((mv = MESH_Next_Vertex(mesh,&idx)))
    MEnt_Rem_AttVal(mv,vidatt);

  idx = 0; i = 0;
  while ((me = MESH_Next_Edge(mesh,&idx)))
    MEnt_Rem_AttVal(me,eidatt);

  idx = 0; i = 0;
  while ((mf = MESH_Next_Face(mesh,&idx)))
    MEnt_Rem_AttVal(mf,fidatt);

  idx = 0; i = 0;
  while ((mr = MESH_Next_Region(mesh,&idx)))
    MEnt_Rem_AttVal(mr,ridatt);
  
  MAttrib_Delete(vidatt);
  MAttrib_Delete(eidatt);
  MAttrib_Delete(fidatt);
  MAttrib_Delete(ridatt);


  fclose(fp);

  return 1;
}

int PMESH_WriteToFile(Mesh_ptr mesh, const char *filename, RepType rtype) {
  FILE *fp;
  char mesg[80], attname[256];
  int i, j, k, idx;
  int gdim, gid;
  int mvid, mvid0, mvid1, mvid2, mrid2, meid, mfid, mrid;
  int nav, nar, nfe, nfv, nrf, nrv, dir=0;
  int nv, ne, nf, nr;
  int natt, ncomp, ival, nent;
  double xyz[3], rval, rdummy, *rval_arr;
  void *pval, *pdummy;
  MVertex_ptr mv, mv0, mv1, mv2;
  MEdge_ptr me;
  MFace_ptr mf;
  MRegion_ptr mr, mr2;
  GEntity_ptr gent;
  List_ptr adjverts, mfedges, mfverts, mrfaces, mrverts, adjregs;
  RepType reptype;
  MAttrib_ptr attrib, vidatt, eidatt, fidatt, ridatt;
  MType attentdim;
  MAttType atttype;
  /*for mpi */
  PType ptype;
  int id, masterparid, globalid;
  /*end for mpi */
  

  if (!(fp = fopen(filename,"w"))) {
    sprintf(mesg,"Cannot open file %-s for writing",filename);
    MSTK_Report("MESH_WriteToFile",mesg,ERROR);
    return 0;
  }

  if (rtype != UNKNOWN_REP) {
    reptype = rtype;
  }
  else {
    reptype = MESH_RepType(mesh);
  }

  nv = MESH_Num_Vertices(mesh);
  ne = MESH_Num_Edges(mesh);
  nf = MESH_Num_Faces(mesh);
  nr = MESH_Num_Regions(mesh);

  fprintf(fp,"MSTK %-2.1lf\n",MSTK_VER);
  fprintf(fp,"%s %d %d %d %d\n",
	  MESH_rtype_str[reptype], 
	  nv, 
	  (reptype >= R1 && reptype <= R4)?0:ne, 
	  (reptype >= R1 && reptype <= R2 && nr)?0:nf, 
	  nr);

  vidatt = MAttrib_New(mesh,"vidatt",INT,MVERTEX);
  eidatt = MAttrib_New(mesh,"eidatt",INT,MEDGE);
  fidatt = MAttrib_New(mesh,"fidatt",INT,MFACE);
  ridatt = MAttrib_New(mesh,"ridatt",INT,MREGION);

  idx = 0; i = 0;
  while ((mv = MESH_Next_Vertex(mesh,&idx)))
    MEnt_Set_AttVal(mv,vidatt,++i,0.0,NULL);

  idx = 0; i = 0;
  while ((me = MESH_Next_Edge(mesh,&idx)))
    MEnt_Set_AttVal(me,eidatt,++i,0.0,NULL);

  idx = 0; i = 0;
  while ((mf = MESH_Next_Face(mesh,&idx)))
    MEnt_Set_AttVal(mf,fidatt,++i,0.0,NULL);

  idx = 0; i = 0;
  while ((mr = MESH_Next_Region(mesh,&idx)))
    MEnt_Set_AttVal(mr,ridatt,++i,0.0,NULL);
  
  fprintf(fp,"vertices\n");
  idx = 0;
  while ((mv = MESH_Next_Vertex(mesh,&idx))) {

    MV_Coords(mv,xyz);

    gdim = MV_GEntDim(mv);
    gid = MV_GEntID(mv);
    /* for mpi */
    
    ptype = MV_PType(mv);
    id = MV_ID(mv);
    masterparid = MV_MasterParID(mv);
    globalid = MV_GlobalID(mv);
    fprintf(fp,"%24.16lf %24.16lf %24.16lf   %d %d\t%12s %6d %3d %6d\n",
	    xyz[0],xyz[1],xyz[2],gdim,gid,MESH_ptype_str[ptype],id,masterparid,globalid);
    /* end for mpi */
    /*    fprintf(fp,"%24.16lf %24.16lf %24.16lf   %d %d\n",
	  xyz[0],xyz[1],xyz[2],gdim,gid); */
    
  }

  if (reptype == R2 || reptype == R4) {
    fprintf(fp,"adjvertices\n");

    idx = 0;
    while ((mv = MESH_Next_Vertex(mesh,&idx))) {

      nav = MV_Num_AdjVertices(mv);
      fprintf(fp,"%d ",nav);
      
      adjverts = MV_AdjVertices(mv);
      for (j = 0; j < nav; j++) {
	mv2 = List_Entry(adjverts,j);
	MEnt_Get_AttVal(mv2,vidatt,&mvid2,&rval,&pval);
	fprintf(fp,"%d ",mvid2);
      }
      fprintf(fp,"\n");
      List_Delete(adjverts);
    }
  }



  if (reptype <= F4 && ne) {
    fprintf(fp,"edges\n");

    idx = 0;
    while ((me = MESH_Next_Edge(mesh,&idx))) {

      mv0 = ME_Vertex(me,0);
      MEnt_Get_AttVal(mv0,vidatt,&mvid0,&rval,&pval);
      mv1 = ME_Vertex(me,1);
      MEnt_Get_AttVal(mv1,vidatt,&mvid1,&rval,&pval);

      gdim = ME_GEntDim(me);
      gid = ME_GEntID(me);
      /* for mpi */
      ptype = ME_PType(me);
      id = ME_ID(me);
      masterparid = ME_MasterParID(me);
      globalid = ME_GlobalID(me);
      fprintf(fp,"%d %d \t%d %d\t%12s %6d %3d %6d\n"
	      ,mvid0,mvid1,gdim,gid,MESH_ptype_str[ptype],id,masterparid,globalid);
      /* end for mpi */
      /* fprintf(fp,"%d %d \t%d %d\n",mvid0,mvid1,gdim,gid); */
    }
  }



  if (reptype <= F4) {

    /* For full representations, always write out faces in terms of edges */

    fprintf(fp,"faces edge\n");
    
    idx = 0;
    while ((mf = MESH_Next_Face(mesh,&idx))) {
      
      nfe = MF_Num_Edges(mf);
      fprintf(fp,"%d ",nfe);
      
      mfedges = MF_Edges(mf,1,0);
      for (j = 0; j < nfe; j++) {
	me = List_Entry(mfedges,j);
	dir = MF_EdgeDir_i(mf,j);
	MEnt_Get_AttVal(me,eidatt,&meid,&rval,&pval);
	if (dir != 1) meid = -meid;
	fprintf(fp,"%d ",meid);
      }
      List_Delete(mfedges);
      
      gdim = MF_GEntDim(mf);
      /*
	gent = MF_GEntity(mf);
	gid = gent ? -99 : 0;
      */
      gid = MF_GEntID(mf);

      /* for mpi */
      ptype = MF_PType(mf);
      id = MF_ID(mf);
      masterparid = MF_MasterParID(mf);
      globalid = MF_GlobalID(mf);
      fprintf(fp,"\t%d %d\t%12s %6d %3d %6d\n",
	      gdim,gid,MESH_ptype_str[ptype],id,masterparid,globalid);
      /* end for mpi */
      /* fprintf(fp,"\t%d %d\n",gdim,gid); */
    }
  }
  else {

    /* For reduced representations, R3 and R4 always write out faces
       in terms of vertices. For reduced representations, R1 and R2
       write out faces in terms of vertices only when there are no
       regions (i.e. faces are the highest level mesh entities) */

    if ((reptype > R2) || (nr == 0)) {

      fprintf(fp,"faces vertex\n");

      idx = 0;
      while ((mf = MESH_Next_Face(mesh,&idx))) {
	
	nfv = MF_Num_Edges(mf);
	fprintf(fp,"%d ",nfv);
	
	mfverts = MF_Vertices(mf,1,0);
	for (j = 0; j < nfv; j++) {
	  mv = List_Entry(mfverts,j);
	  MEnt_Get_AttVal(mv,vidatt,&mvid,&rval,&pval);
	  fprintf(fp,"%d ",mvid);
	}
	List_Delete(mfverts);

	gdim = MF_GEntDim(mf);
	gid = MF_GEntID(mf);
	/* for mpi */
	ptype = MF_PType(mf);
	id = MF_ID(mf);
	masterparid = MF_MasterParID(mf);
	globalid = MF_GlobalID(mf);
	fprintf(fp,"\t%d %d\t%12s %6d %3d %6d\n",
		gdim,gid,MESH_ptype_str[ptype],id,masterparid,globalid);
	/* end for mpi */
	/* fprintf(fp,"\t%d %d\n",gdim,gid);*/
      }
    }
	
  }


  if (nr) {
    if (reptype <= F4 || reptype >= R2) {
      fprintf(fp,"regions face\n");

      idx = 0;
      while ((mr = MESH_Next_Region(mesh,&idx))) {

	nrf = MR_Num_Faces(mr);
	fprintf(fp,"%d ",nrf);

	mrfaces = MR_Faces(mr);
	for (j = 0; j < nrf; j++) {
	  mf = List_Entry(mrfaces,j);
	  dir = MR_FaceDir_i(mr,j);
	  MEnt_Get_AttVal(mf,fidatt,&mfid,&rval,&pval);
	  if (dir != 1) mfid = -mfid;
	  fprintf(fp,"%d ",mfid);
	}
	List_Delete(mrfaces);
	
	gdim = MF_GEntDim(mr);
	gid = MR_GEntID(mr);
	/* for mpi */
	ptype = MR_PType(mr);
	id = MR_ID(mr);
	masterparid = MR_MasterParID(mr);
	globalid = MR_GlobalID(mr);
	fprintf(fp,"\t%d %d\t%12s %6d %3d %6d\n",
		gdim,gid,MESH_ptype_str[ptype],id,masterparid,globalid);
	/* end for mpi */
	/* fprintf(fp,"\t%d %d\n",gdim,gid); */
      }
    }
    else {
      fprintf(fp,"regions vertex\n");

      idx = 0;
      while ((mr = MESH_Next_Region(mesh,&idx))) {

	nrv = MR_Num_Vertices(mr);
	fprintf(fp,"%d ",nrv);

	mrverts = MR_Vertices(mr);
	for (j = 0; j < nrv; j++) {
	  mv = List_Entry(mrverts,j);
	  MEnt_Get_AttVal(mv,vidatt,&mvid,&rval,&pval);
	  fprintf(fp,"%d ",mvid);
	}
	List_Delete(mrverts);
	
	gdim = MR_GEntDim(mr);
	gid = MR_GEntID(mr);
	/* for mpi */
	ptype = MR_PType(mr);
	id = MR_ID(mr);
	masterparid = MR_MasterParID(mr);
	globalid = MR_GlobalID(mr);
	fprintf(fp,"\t%d %d\t%12s %6d %3d %6d\n",
		gdim,gid,MESH_ptype_str[ptype],id,masterparid,globalid);
	/* end for mpi */

	/*fprintf(fp,"\t%d %d\n",gdim,gid);*/
      }
    }

    if (reptype == R2 || reptype == R4) {
      fprintf(fp,"adjregions\n");
      
      idx = 0;
      while ((mr = MESH_Next_Region(mesh,&idx))) {

	nar = MR_Num_Faces(mr);
	fprintf(fp,"%d ",nar);

	adjregs = MR_AdjRegions(mr);

	for (j = 0; j < nar; j++) {
	  mr2 = List_Entry(adjregs,j);
	  if ((long) mr2 == -1) 
	    fprintf(fp,"%d ",0);
	  else {
	    MEnt_Get_AttVal(mr2,ridatt,&mrid2,&rval,&pval);
	    fprintf(fp,"%d ",mrid2);
	  }
	}
	fprintf(fp,"\n");
	List_Delete(adjregs);
      }
    }
  }


  /* Write out attributes if there are more than the 4 that we created 
    in this routine */


  if ((natt = MESH_Num_Attribs(mesh)) > 4) {

    fprintf(fp,"attributes\n");

    for (i = 0; i < natt; i++) {
      
      attrib = MESH_Attrib(mesh,i);

      /* Don't write out attribs we created for the internal use of 
	 this routine */
      if (attrib == vidatt || attrib == eidatt || attrib == fidatt || 
	  attrib == ridatt) continue;
      
      MAttrib_Get_Name(attrib,attname);

      atttype = MAttrib_Get_Type(attrib);
      if (atttype == POINTER) continue;  /* cannot write it out */

      ncomp = MAttrib_Get_NumComps(attrib);

      attentdim = MAttrib_Get_EntDim(attrib);


      /* First count how many entities actually have the attribute assigned */

      nent = 0;
      switch(attentdim) {
      case MVERTEX:
	idx = 0;
	while ((mv = MESH_Next_Vertex(mesh,&idx)))
	  if (MEnt_Get_AttVal(mv,attrib,&ival,&rval,&pval)) nent++;
	break;
      case MEDGE:
	idx = 0;
	while ((me = MESH_Next_Edge(mesh,&idx)))
	  if (MEnt_Get_AttVal(me,attrib,&ival,&rval,&pval)) nent++;
	break;
      case MFACE:
	idx = 0;
	while ((mf = MESH_Next_Face(mesh,&idx)))
	  if (MEnt_Get_AttVal(mf,attrib,&ival,&rval,&pval)) nent++;	    
	break;
      case MREGION: 
	idx = 0;
	while ((mr = MESH_Next_Region(mesh,&idx)))
	  if (MEnt_Get_AttVal(mr,attrib,&ival,&rval,&pval)) nent++;
	break;
      case MALLTYPE:
	idx = 0;
	while ((mv = MESH_Next_Vertex(mesh,&idx)))
	  if (MEnt_Get_AttVal(mv,attrib,&ival,&rval,&pval)) nent++;
	idx = 0;
	while ((me = MESH_Next_Edge(mesh,&idx)))
	  if (MEnt_Get_AttVal(me,attrib,&ival,&rval,&pval)) nent++;
	idx = 0;
	while ((mf = MESH_Next_Face(mesh,&idx)))
	  if (MEnt_Get_AttVal(mf,attrib,&ival,&rval,&pval)) nent++;	    
	idx = 0;
	while ((mr = MESH_Next_Region(mesh,&idx)))
	  if (MEnt_Get_AttVal(mr,attrib,&ival,&rval,&pval)) nent++;
	break;	
      default:
	break;
      } /* switch (attentdim) */


      /* No point in writing out attribute if no entity uses it! Or is there? */

      if (!nent) continue;



      fprintf(fp,"%-s\n",attname);

      switch(atttype) {
      case INT:
	fprintf(fp,"INT\n");
	break;
      case DOUBLE:
	fprintf(fp,"DOUBLE\n");
	break;
      case VECTOR:
	fprintf(fp,"VECTOR\n");
	break;
      case TENSOR:
	fprintf(fp,"TENSOR\n");
	break;
      default:
	MSTK_Report("MESH_WriteToFile",
		    "Unrecognizable or unprintable attribute type\n",WARN);
	continue;	
      }

      fprintf(fp,"%-d\n",ncomp);

      switch(attentdim) {
      case MVERTEX:
	fprintf(fp,"MVERTEX\n");
	break;
      case MEDGE:
	fprintf(fp,"MEDGE\n");
	break;
      case MFACE:
	fprintf(fp,"MFACE\n");
	break;
      case MREGION:
	fprintf(fp,"MREGION\n");
	break;
      case MALLTYPE:
	fprintf(fp,"MALLTYPE\n");
	break;
      default:
	MSTK_Report("Mesh_WriteToFile","Unrecognized entity type",WARN);
	break;
      }

      fprintf(fp,"%-d\n",nent);


      switch(attentdim) {
      case MVERTEX:
	idx = 0;
	while ((mv = MESH_Next_Vertex(mesh,&idx))) {
	  if (MEnt_Get_AttVal(mv,attrib,&ival,&rval,&pval)) {
	    MEnt_Get_AttVal(mv,vidatt,&mvid,&rdummy,&pdummy);
	    fprintf(fp,"0 %-d ",mvid);
	    switch (atttype) {
	    case INT:
	      fprintf(fp," %-d",ival);
	      break;
	    case DOUBLE: 
	      fprintf(fp," %-lf ",rval);
	      break;
	    case VECTOR: case TENSOR:
	      rval_arr = (double *) pval;
	      for (k = 0; k < ncomp; k++)
		fprintf(fp," %-lf ",rval_arr[k]);
	      break;
	    default:
	      break;
	    }
	    fprintf(fp,"\n");
	  }
	}
	break;
      case MEDGE:
	idx = 0;
	while ((me = MESH_Next_Edge(mesh,&idx))) {
	  if (MEnt_Get_AttVal(me,attrib,&ival,&rval,&pval)) {
	    MEnt_Get_AttVal(me,eidatt,&meid,&rdummy,&pdummy);
	    fprintf(fp,"1 %-d ",meid);
	    switch (atttype) {
	    case INT:
	      fprintf(fp," %-d",ival);
	      break;
	    case DOUBLE: 
	      fprintf(fp," %-lf ",rval);
	      break;
	    case VECTOR: case TENSOR:
	      rval_arr = (double *) pval;
	      for (k = 0; k < ncomp; k++)
		fprintf(fp," %-lf ",rval_arr[k]);
	      break;
	    default:
	      break;
	    }
	    fprintf(fp,"\n");
	  }
	}
	break;
      case MFACE:
	idx = 0;
	while ((mf = MESH_Next_Face(mesh,&idx))) {
	  if (MEnt_Get_AttVal(mf,attrib,&ival,&rval,&pval)) {
	    MEnt_Get_AttVal(mf,fidatt,&mfid,&rdummy,&pdummy);
	    fprintf(fp,"2 %-d ",mfid);
	    switch (atttype) {
	    case INT:
	      fprintf(fp," %-d",ival);
	      break;
	    case DOUBLE: 
	      fprintf(fp," %-lf ",rval);
	      break;
	    case VECTOR: case TENSOR:
	      rval_arr = (double *) pval;
	      for (k = 0; k < ncomp; k++)
		fprintf(fp," %-lf ",rval_arr[k]);
	      break;
	    default:
	      break;
	    }
	    fprintf(fp,"\n");
	  }
	}
	break;

      case MREGION: 
	idx = 0;
	while ((mr = MESH_Next_Region(mesh,&idx))) {
	  if (MEnt_Get_AttVal(mr,attrib,&ival,&rval,&pval)) {
	    MEnt_Get_AttVal(mr,ridatt,&mrid,&rdummy,&pdummy);
	    fprintf(fp,"3 %-d ",mrid);
	    switch (atttype) {
	    case INT:
	      fprintf(fp," %-d",ival);
	      break;
	    case DOUBLE: 
	      fprintf(fp," %-lf ",rval);
	      break;
	    case VECTOR: case TENSOR:
	      rval_arr = (double *) pval;
	      for (k = 0; k < ncomp; k++)
		fprintf(fp," %-lf ",rval_arr[k]);
	      break;
	    default:
	      break;
	    }
	    fprintf(fp,"\n");
	  }
	}
	break;

      case MALLTYPE:
	idx = 0;
	while ((mv = MESH_Next_Vertex(mesh,&idx))) {
	  if (MEnt_Get_AttVal(mv,attrib,&ival,&rval,&pval)) {
	    MEnt_Get_AttVal(mv,vidatt,&mvid,&rdummy,&pdummy);
	    fprintf(fp,"0 %-d ",mvid);
	    switch (atttype) {
	    case INT:
	      fprintf(fp," %-d",ival);
	      break;
	    case DOUBLE: 
	      fprintf(fp," %-lf ",rval);
	      break;
	    case VECTOR: case TENSOR:
	      rval_arr = (double *) pval;
	      for (k = 0; k < ncomp; k++)
		fprintf(fp," %-lf ",rval_arr[k]);
	      break;
	    default:
	      break;
	    }
	    fprintf(fp,"\n");
	  }
	}
	idx = 0;
	while ((me = MESH_Next_Edge(mesh,&idx))) {
	  if (MEnt_Get_AttVal(me,attrib,&ival,&rval,&pval)) {
	    MEnt_Get_AttVal(me,eidatt,&meid,&rdummy,&pdummy);
	    fprintf(fp,"1 %-d ",meid);
	    switch (atttype) {
	    case INT:
	      fprintf(fp," %-d",ival);
	      break;
	    case DOUBLE: 
	      fprintf(fp," %-lf ",rval);
	      break;
	    case VECTOR: case TENSOR:
	      rval_arr = (double *) pval;
	      for (k = 0; k < ncomp; k++)
		fprintf(fp," %-lf ",rval_arr[k]);
	      break;
	    default:
	      break;
	    }
	    fprintf(fp,"\n");
	  }
	}
	idx = 0;
	while ((mf = MESH_Next_Face(mesh,&idx))) {
	  if (MEnt_Get_AttVal(mf,attrib,&ival,&rval,&pval)) {
	    MEnt_Get_AttVal(mf,fidatt,&mfid,&rdummy,&pdummy);
	    fprintf(fp,"2 %-d ",mfid);
	    switch (atttype) {
	    case INT:
	      fprintf(fp," %-d",ival);
	      break;
	    case DOUBLE: 
	      fprintf(fp," %-lf ",rval);
	      break;
	    case VECTOR: case TENSOR:
	      rval_arr = (double *) pval;
	      for (k = 0; k < ncomp; k++)
		fprintf(fp," %-lf ",rval_arr[k]);
	      break;
	    default:
	      break;
	    }
	    fprintf(fp,"\n");
	  }
	}
	idx = 0;
	while ((mr = MESH_Next_Region(mesh,&idx))) {
	  if (MEnt_Get_AttVal(mr,attrib,&ival,&rval,&pval)) {
	    MEnt_Get_AttVal(mr,ridatt,&mrid,&rdummy,&pdummy);
	    fprintf(fp,"3 %-d ",mrid);
	    switch (atttype) {
	    case INT:
	      fprintf(fp," %-d",ival);
	      break;
	    case DOUBLE: 
	      fprintf(fp," %-lf ",rval);
	      break;
	    case VECTOR: case TENSOR:
	      rval_arr = (double *) pval;
	      for (k = 0; k < ncomp; k++)
		fprintf(fp," %-lf ",rval_arr[k]);
	      break;
	    default:
	      break;
	    }
	    fprintf(fp,"\n");
	  }
	}
	break;	
      default:
	break;
      } /* switch (attentdim) */

    } /* for (i = 0; i < natt) */
    
  } /* if (Mesh_Num_Attribs(mesh)) */
  

  idx = 0; i = 0;
  while ((mv = MESH_Next_Vertex(mesh,&idx)))
    MEnt_Rem_AttVal(mv,vidatt);

  idx = 0; i = 0;
  while ((me = MESH_Next_Edge(mesh,&idx)))
    MEnt_Rem_AttVal(me,eidatt);

  idx = 0; i = 0;
  while ((mf = MESH_Next_Face(mesh,&idx)))
    MEnt_Rem_AttVal(mf,fidatt);

  idx = 0; i = 0;
  while ((mr = MESH_Next_Region(mesh,&idx)))
    MEnt_Rem_AttVal(mr,ridatt);
  
  MAttrib_Delete(vidatt);
  MAttrib_Delete(eidatt);
  MAttrib_Delete(fidatt);
  MAttrib_Delete(ridatt);


  fclose(fp);

  return 1;
}




#ifdef __cplusplus
}
#endif
